FFT For Image Analysis

Gati Aher
SP2021: Quantitative Engineering Analysis

Explore Image Fourier Transforms

sine plots
explore_fft(imread("img/vertical_sine_img_f8.png"));
explore_fft(imread("img/horizontal_sine_img_f8.png"));
explore_fft(imread("img/diagonal_sine_img_scaled_f8.png"));
explore_fft(imread("img/vertical_sine_img_f32.png"));
explore_fft(imread("img/horizontal_sine_img_f32.png"));
explore_fft(imread("img/diagonal_sine_img_scaled_f32.png"));
rectangles
explore_fft(imread("img/rect_32x64px.png"));
explore_fft(imread("img/rect_32x64px_offset_by_50.png"));
explore_fft(imread("img/square_32px.png"));
explore_fft(imread("img/square_64px.png"));
explore_fft(imread("img/square_128px.png"));
natural patterns
explore_fft(imread("img/brick.png"));
explore_fft(imread("img/honeycomb.png"));
explore_fft(imread("img/mustard_seeds.png"));
complex images
explore_fft(imread("img/cameraman.png"));
explore_fft(imread("img/sugar_cube.png"));
explore_fft(imread("img/flower.png"));
explore_fft(imread("img/bacteria_slayer.png"));
other
explore_fft("img/blank.png");
explore_fft("img/diagonal_lines.png");
explore_fft("img/noise.png");

Windowing

img = imread("img/diagonal_sine_img_scaled_f32.png");
w = hann(N, 'periodic');
explore_fft(double(img) .* (w.*w'));
img = imread("img/flower.png");
w = hann(N, 'periodic');
explore_fft(double(img) .* (w.*w'));
img = imread("img/bacteria_slayer.png");
w = hann(N, 'periodic') .* 0.25;
explore_fft(double(img) .* (w.*w'));

Check Filters

Sanity check that filters work as expected:
N = 256;
figure, imshow(distance_mat(zeros(N, N)), [0, 100]), title("Distance Matrix");
idf = idealFilterLP(N/4, zeros(N, N));
bwf5 = butterworthLP(N/4, zeros(N, N), 5);
bwf2 = butterworthLP(N/4, zeros(N, N), 2);
gf = gaussianLP(N/4, zeros(N, N));
subplot(2, 4, 1), imshow(idf, []), title("Ideal Low Pass Filter");
subplot(2, 4, 2), imshow(bwf5, []), title("Butterworth Order 5 Low Pass Filter");
subplot(2, 4, 3), imshow(bwf2, []), title("Butterworth Order 2 Low Pass Filter");
subplot(2, 4, 4), imshow(gf, []), title("Gaussian Low Pass Filter");
subplot(2, 4, 5), imshow(idealFilterHP(N/4, zeros(N, N)), [0, 1]), title("Ideal High Pass Filter");
subplot(2, 4, 6), imshow(butterworthHP(N/4, zeros(N, N), 5), [0, 1]), title("Butterworth Order 5 High Pass Filter");
subplot(2, 4, 7), imshow(butterworthHP(N/4, zeros(N, N), 2), [0, 1]), title("Butterworth Order 2 High Pass Filter");
subplot(2, 4, 8), imshow(gaussianHP(N/4, zeros(N, N)), [0, 1]), title("Gaussian High Pass Filter");
% calculate fft phase and magnitude for each filter
explore_fft(idf);
explore_fft(bwf5);
explore_fft(bwf2);
explore_fft(gf);

Try Filtering

img = imread("img/flower.png");
N = size(img, 1);
explore_fft(img);
explore_filter_fft(img, idealFilterLP(N/4, zeros(N, N)));
Warning: Displaying real part of complex input.
explore_filter_fft(img, butterworthLP(N/4, zeros(N, N), 5));
Warning: Displaying real part of complex input.
explore_filter_fft(img, butterworthLP(N/4, zeros(N, N), 2));
Warning: Displaying real part of complex input.
explore_filter_fft(img, gaussianLP(N/4, zeros(N, N)));
Warning: Displaying real part of complex input.
figure
subplot(2, 5, 1), imshow(img), title("Original Image");
subplot(2, 5, 2), imshow(get_filter_fft(img, idealFilterHP(N/4, zeros(N, N))), []), title("Recon. Ideal HP size N/4");
Warning: Displaying real part of complex input.
subplot(2, 5, 3), imshow(get_filter_fft(img, butterworthHP(N/4, zeros(N, N), 5)), []), title("Recon. Butterworth 5 HP size N/4");
Warning: Displaying real part of complex input.
subplot(2, 5, 4), imshow(get_filter_fft(img, butterworthHP(N/4, zeros(N, N), 2)), []), title("Recon. Butterworth 2 HP size N/4");
Warning: Displaying real part of complex input.
subplot(2, 5, 5), imshow(get_filter_fft(img, gaussianHP(N/4, zeros(N, N))), []), title("Recon. Gaussian HP size N/4");
Warning: Displaying real part of complex input.
subplot(2, 5, 7), imshow(get_filter_fft(img, idealFilterHP(N/32, zeros(N, N))), []), title("Recon. Ideal HP size N/32");
Warning: Displaying real part of complex input.
subplot(2, 5, 8), imshow(get_filter_fft(img, butterworthHP(N/32, zeros(N, N), 5)), []), title("Recon. Butterworth 5 HP size N/32");
Warning: Displaying real part of complex input.
subplot(2, 5, 9), imshow(get_filter_fft(img, butterworthHP(N/32, zeros(N, N), 2)), []), title("Recon. Butterworth 2 HP size N/32");
Warning: Displaying real part of complex input.
subplot(2, 5, 10), imshow(get_filter_fft(img, gaussianHP(N/32, zeros(N, N))), []), title("Recon. Gaussian HP size N/32");
Warning: Displaying real part of complex input.

FFT Function:

function f = explore_fft(img)
img_fft = fft2(double(img));
img_phase = angle(img_fft);
img_mag = abs(img_fft);
img_mag_log = log(1 + img_mag);
img_fftshift = fftshift(img_fft);
img_phase_shift = angle(img_fftshift);
img_mag_shift = abs(img_fftshift);
img_mag_shift_log = log(1 + img_mag_shift);
f = figure;
subplot(2, 4, 1), imshow(img), title("Original Image");
subplot(2, 4, 2), imshow(img_phase, []), title("Phase Plot");
subplot(2, 4, 3), imshow(img_mag, []), title("Magnitude Plot");
subplot(2, 4, 4), imshow(img_mag_log, []), title("Log of Magnitude");
subplot(2, 4, 6), imshow(img_phase_shift, []), title("Shifted Phase Plot");
subplot(2, 4, 7), imshow(img_mag_shift, []), title("Shifted Magnitude Plot");
subplot(2, 4, 8), imshow(img_mag_shift_log, []), title("Log of Shifted Magnitude");
end

FFT + Filter Function

function f = explore_filter_fft(img, filter)
img_fftshift = fftshift(fft2(double(img)));
img_filter_fft = img_fftshift .* double(filter);
img_filter_decentralize = ifftshift(img_filter_fft);
img_filter_reconstructed = ifft2(img_filter_decentralize);
f = figure;
subplot(2, 3, 1), imshow(filter), title("Filter");
subplot(2, 3, 2), imshow(img), title("Original Image");
subplot(2, 3, 3), imshow(log(1 + abs(img_fftshift)), []), title("Log of Centralized Magnitude Plot");
subplot(2, 3, 4), imshow(log(1 + abs(img_filter_fft)), []), title("Log of Filtered Centralized Magnitude Plot");
subplot(2, 3, 5), imshow(log(1 + abs(img_filter_decentralize)), []), title("Log of Decentralized Magnitude Plot");
subplot(2, 3, 6), imshow(img_filter_reconstructed, []), title("Reconstructed Image");
end
function img_filter_reconstructed = get_filter_fft(img, filter)
img_fftshift = fftshift(fft2(double(img)));
img_filter_fft = img_fftshift .* double(filter);
img_filter_decentralize = ifftshift(img_filter_fft);
img_filter_reconstructed = ifft2(img_filter_decentralize);
end

Filter Functions:

function distMat = distance_mat(img)
% get center of image
[nx, ny] = size(img);
x = linspace(-nx/2, nx/2, nx);
y = linspace(-ny/2, ny/2, ny);
% calculate distance matrix
[X, Y] = meshgrid(x, y);
distMat = sqrt((X.^2)+(Y.^2));
end
function base = idealFilterLP(cutoff, img)
% get distance matrix for image
distMat = distance_mat(img);
% perform logical mask
mask = distMat < cutoff;
% everything below cutoff is True (1)
% everything above cutoff is False (0)
base = mask;
end
function base = idealFilterHP(cutoff, img)
% get distance matrix for image
distMat = distance_mat(img);
% perform logical mask
mask = distMat > cutoff;
% everything below cutoff is False (0)
% everything above cutoff is True (1)
base = mask;
end
function base = butterworthLP(cutoff, img, order)
% get distance matrix for image
distMat = distance_mat(img);
% get butterworth values
base = 1./(1+(distMat./cutoff).^(2*order));
end
function base = butterworthHP(cutoff, img, order)
% high-pass is just inverted low-pass
base = 1 - butterworthLP(cutoff, img, order);
end
function base = gaussianLP(cutoff, img)
% get distance matrix for image
distMat = distance_mat(img);
% get butterworth values
% base = expm(((-1 .* (distMat.^2))./(2*(cutoff^2))));
base = exp(-1 .* (distMat.^2)./(2*(cutoff^2)));
end
function base = gaussianHP(cutoff, img)
% high-pass is just inverted low-pass
base = 1 - gaussianLP(cutoff, img);
end